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Hybrid Molecular Dynamics simulations of living filaments 

Mathieu Caby/ Priscilla Hardas/ Sanoop Ramachandran,-'^ and Jean-Paul Rvckaert-'^'Rl 

Physique des Polymeres, Universite Libre de Bruxelles, 
Campus Plaine, CP 223, B-1050 Brussels, Belgium 

(Dated: 29 February 2012) 

We propose a hybrid Molecular Dynamics/Multi-particle Collision Dynamics model to simulate a set of 
self-assembled semiflexible filaments and free monomers. Further, we introduce a Monte-Carlo scheme to 
deal with single monomer addition (polymerization) or removal (depolymerization) , satisfying the detailed 
balance condition within a proper statistical mechanical framework. This model of filaments, based on the 
wormlike chain, aims to represent equilibrium polymers with distinct reaction rates at both ends, such as 
self-assembled ADP-actin filaments in the absence of ATP hydrolysis and other proteins. We report the 
distribution of filament lengths and the corresponding dynamical fluctuations on an equilibrium trajectory. 
Potential generalizations of this method to include irreversible steps like ATP-actin hydrolysis are discussed. 



I. INTRODUCTION 

"Living" polymers form a class of supramolecular 
structures characterized by the reversible self-assembly 
of monomeric units into chains at equilibriumiir— Some 
examples are systems forming discotic structures^i^ fib- 
rillar structures comprising of oligo(p-phenylenevinylene) 
derivatives)^ and wormlike micellesi^ Filamentous struc- 
tures of proteins found in the cell, such as F-actin^SiiiS 
microtubules,—'^ intermediate filaments^^ and bacte- 
rial piliiiii^ are usually coined as non-equilibrium poly- 
mers Their assembly /disassembly properties are cou- 
pled to irreversible chemical steps, like the ATP/GTP 
(adenosine triphosphate/guanosine triphosphate) hydrol- 
ysis to ADP/GDP (adenosine diphosphate/guanosine 
diphosphate) within the filament complexes which act 
as catalysts. These complexes containing hydrolyzed nu- 
cleotides turn out to be more loosely bounded within 
the filaments. This leads to a shortening of persis- 
tence length for the hydrolyzed portions and to differ- 
ent critical concentrations at active ends. In contrast 
to the passive polymerization of equilibrium polymers, 
the active polymerization/dcpolymerization confers to 
these polar biofilaments unusual properties such as tread- 
milling where the filament grows at one end and shrinks 
at the other, or such as generation of work at the ex- 
pense of chemical energy when a temporary network of 
filaments pushes on a membrane (e.g. in filopodia),^'^ 
Numerous in vitro experiments have been conducted to 
observe and analyze various peculiar properties of biofil- 
aments. Tubulin fibers have been observed to alternate 
between slow growth and catastrophe shrinkage — al- 
lowing for rapid reorganization of the cytoskeleton net- 
work while the role played by geometric boundaries on 
the growth dynamics of actin bundles was highlighted 
recently by Reymann et alJ^ 

There have been several theoretical investigations 
based on mean-field stochastic models to study the con- 
tour length and composition dynamics of single free actin 



(also microtubule) filaments in presence of ATP hydroly- 
These theories have been adapted to investigate 
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the behavior of bundles of filaments pushing against a 
flexible membrane^^ or wall^i using straight and rigid 
filaments. They treat the free monomer concentration 
as homogeneous in space and constant in time and sim- 
plified hypotheses are assumed on the load distribution 
exerted by the surface on the filaments. Under condi- 
tions of confinement or crowding of bundles, it is expected 
that spatial correlations between filaments, their flexibil- 
ity, the free monomer concentration fluctuations and the 
presence of local flows interfere with the main polymer- 
ization/depolymerization mechanism. To deal with more 
sophisticated situations where these aspects could be si- 
multaneously taken into account, mesoscopic simulation 
approaches are needed. Pioneering Brownian Dynamics 
(BD) simulations, treating at the same time monomer 
diffusion and filament/network self-assembly, have been 
reported recently. They simulate actin networks push- 
ing against an objec t^^i^^ or a long single actin filament 
undergoing (de)polymerization steps in presence of ATP 
hydrolysis j^2i2^ 

In this paper, we propose and develop a mesoscopic 
model to simulate a set of self-assembled semiflexi- 
ble filaments within a statistical mechanical framework. 
The total number of fllaments and the total number of 
monomers are both fixed in the simulations, leading to 
semi-grand canonical ensemble at equilibrium. The fila- 
ments are modeled according to a discretized version of 
the wormlike-chain with a variable contour length and 
an adjustable persistence length. The chemical reactive 
steps are simulated via a novel Monte-Carlo (MC) move 
which forms the core of the present work. This MC 
move which satisfies the micro-reversibility conditions, 
is specifically suitable to model single addition/removal 
of a monomer at the active ends of the semiflexible flla- 
mcnt. It is reminiscent of the MC moves developed for 
modeling scission/recombination kinetics in flexible lin- 
ear micelles^i^S or to model networks of reversible associ- 
ating polymers^ in which operational parameters can be 
fixed to separately adjust the equilibrium constant and 
the barrier height for any reversible chemical reaction. In 
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FIG. 1. A scliematic showing tlie model filament with beads 
of size G and bond length Iq. The angle between successive 
bonds are denoted by 9. 



absence of any additional irreversible step, we get a model 
for a set of equilibrium scmiflcxiblc filaments in a well de- 
fined ensemble. All these modeling aspects are detailed 
in Section ini Section HITl covers the theoretical treatment 
of this system, starting with an ideal mixture treatment 
and looking then for packing effects on the filament size 
distribution and on the kinetic rates. Our simulation ex- 
periments are detailed and their results analyzed in Sec- 
tions HV] and |V] respectively. In the Section |Vl], we come 
back on some methodological issues and envisage the ex- 
tension to non-equilibrium polymers which requires the 
consideration of monomer flags to distinguish between 
hydrolyzed and non-hydrolyzed ATP complexes and the 
consideration of additional chemical steps. The applica- 
tion of our methodology to a specific biofilament is dis- 
cussed by providing the explicit parametrization which 
would be needed to model ADP-F-actin (homopolymer of 
ADP-actin complexes in interaction with free ADP-actin 
monomers). This is an interesting example, not found 
in vivo, of a polar equilibrium polymer where the poly- 
merization and depolymerization steps take place at two 
kinetically incqual ends characterized by the same critical 
concentration. The paper closes with a short conclusion 
in Section IVlIl 



A. The solute 

The solute system treated throughout this work is a 
set of Nt spherical monomers of mass m. A fraction 
of these are present as free monomers. The remaining 
monomers constitute the building blocks oi N f assembled 
polydisperse filaments with a range of sizes going from 
a minimum of three monomers to a maximum (set to 
an integer value z). The justification and the general 
consequences of imposed boundaries to the filament size 
will be discussed in Sections III CI and IIIDI 

These filaments are modeled as semiflexible polymers 
with excluded volume interactions with its basic fea- 
tures adapted from a standard modeli^"— A schematic 
of the filament is shown in Fig. [TJ Any filament with i 
monomers has i — 1 bonding interactions Ui{r) with an 
equilibrium length Iq and an energy depth of — eo 



1 . 



U,{T) = -k{r-lof -e, = Ul'^ 



(r) + Ul 



Icc 



(1) 



in terms of the absolute distance r between adjacent 
monomers. The first purely vibrational term UJ^^{r), 
where k is the stretching force constant, maintains a 
chain structure with contour length Lc — {i — l)^o- The 
second "electronic" term — eo corresponds to a (negative) 
bonding energy with respect to a reference zero energy 
level corresponding to the large distance non-covalent 
pair energy between monomers, in order to favor the 
polymerization step (sec Section III C|) . 

These non-covalent interactions are introduced via a 
shifted and truncated Lennard-Jones (LJ) potential, de- 
noted by U2{r). 



U2{r) = 4e 



e(2i/V-r). (2) 



Here e and a are the energy and diameter parameters of 
the LJ potential. The Heaviside step function satisfies 



e(x) 



for x < 0, 

1 for a; > 0. 



(3) 



II. GENERAL MESOSCOPIC MODEL OF A SET OF 
INTERACTING LIVING FILAMENTS AND FREE 
MONOMERS IN A THERMAL BATH 



We first describe the detailed microscopic model of 
the solute/solvent bath system. Our simulation method 
consists of four parts: the adopted solute model for a 
mixture of free monomers and semiflexible filaments of 
various lengths (Section III Ap . the solvent bath in which 
the solutes are immersed and the nature of the solute- 
solvent coupling (Section IIIB[) . an original stochastic 
(de)polymerization procedure which, together with the 
solute model, defines the "living" filaments concept (Sec- 
tion lll C|) and the statistical mechanics ensemble in which 
our finite-size system is sampled by our hybrid simula- 
tions (Section IIIPp . 



Both intramolecular (beyond second neighbors) and in- 
termolecular pairs are subject to such excluded volume 
interactions. Finally, a set of i — 2 three-body poten- 
tials (j)^'^'^'^ {0) is used to favor the straight alignment of 
successive bonds [6 — 0) within each filament 



/bend 



[e) = /s'(i -cos( 



(4) 



where k' is the bending stiffness energy parameter, di- 
rectly related to the persistence length Zp according to 

k' = kBTlp/lo. 

Note that in applications of the model to filaments, one 
typically has Ip 3> and therefore, the short range (in 
physical distance) of the repulsive non-covalent interac- 
tions U2{r) implies that the latter will usually not operate 
as "local interactions" but possibly as "long-range inter- 
actions" if the filaments are sufficiently long to behave as 
coils(Lc ^ Ip). 



Simulations of living filaments 

The particles trajectories are obtained by solving the 
Newton's equations using the velocity- Verlet integrator, 
as in any standard Molecular Dynamics (MD) simulation. 

B. The solvent bath 

The system of Nt monomers is immersed in a sea 
of solvent particles of mass treated according to 
the Multi-Particle CoUisional Dynamics (MPCD) model. 
This technique, first introduced by Malevancts and 
Kapral^i^^ uses simplified solvent dynamics that faith- 
fully reproduces the correct long time fiuid behavior. 
Subsequent to its introduction, it has been extensively 
studied theoretically as well as being used to investigate 
various complex systems such as colloids, polymers and 
vesicles immersed in a solvent Our use of this solute- 
solvent model, except for the living filaments aspects, is 
close to the semiflexible chains in solution.— 

In this solvent model, all solvent particles move as free 
particles between (local) collisions taking place every At 
time steps. The simulation volume V is divided into cubic 
cells of side uq. The local collisions which take place 
independently in each cell imply velocity changes of all 
particles of the cell (solvent particles, free monomers and 
monomers pertaining to a filament) while preserving the 
total mass, linear momentum and energy of that cell. 
Explicitly, within each cell, the relative velocity of each 
particle (with respect to the cell subsystem center of mass 
velocity) is rotated by a fixed angle a around a direction 
selected at random on a sphere with uniform probability. 

Between the collision steps, the solute dynamics is fol- 
lowed by solving the equations of motion of the solute 
subsystem, treated independently of the solvent. This 
coupling enables the proper thermalization of the solute 
particles as well as the incorporation of hydrodynamic 
effects. 



C. Modeling of (de)polymerlzatlon steps 

A fundamental property of bio-filamcnts such as actin, 
tubulin, pili etc., is their ability to grow/shrink via poly- 
merization/depolymerization reactions at their ends.— 
We propose to model these chemical reactions as in- 
stantaneous events which modify locally the topological 
connections of the so-called "reacting monomer". This 
change will be governed by a MC scheme, satisfying de- 
tailed balance, which is fully consistent with all the topo- 
logical arrangements allowed by the semi-grand ensem- 
ble equilibrium partition function of monomers/filaments 
reacting mixture. Any instantaneous chemical reaction 
event taking place in our system will connect two micro- 
scopic states / and J differing only by the position of the 
single reacting monomer of (say index ii) and by topo- 
logical and geometrical constraints involving its link with 
the last monomer, (say of index 12) of a particular active 
end of a filament and its immediate neighbor in the fil- 
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FIG. 2. State I: The polymerized state where the react- 
ing monomer shown in dark blue is connected to the fila- 
ment through intramolecular interactions and is, necessarily, 
located in a low intramolecular energy volume resulting from 
the intersection of a conic volume (in green) and the spher- 
ical layer (red) which respectively limit the bending and the 
stretching energies (see main text). The dashed blue circles 
represent the volume Vs and the light green particles represent 
free monomers. 




FIG. 3. State J: The depolymerized state where the react- 
ing monomer shown in blue is free but is necessarily iocated 
within the spherical layer shown with daslied-blue lines wliose 
minimum radius limits the excluded volume pair interaction 
energy with the fast monomer of the target filament end (see 
main text). 

ament sequence (say index 13). These three monomers 
which belong to the global set of Nt solute monomers, 
are directly implied in the chemical step which we now 
describe with the help of Figs. [5] and [3] respectively rep- 
resenting the initial and final configurations of a poly- 
merized or depolymerized state. In these figures, the 
dashed bead indicates the position of the unique reacting 
monomer ii in the alternative state. 

In the polymerized state / the reacting monomer, 
shaded dark blue in Fig. [51 is linked to the last monomer 
i2 of the filament through a Ui{r) bonding potential and 
contributes to the bending energy term (j)^'''^'^{9) with 
monomers 12 and 13. It is further imposed that in state /, 
the polymerized monomer of index ii strictly lies within 
the region of low intramolecular energy volume Vc deter- 
mined as the intersection of two volumes: i) the conical 
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volume whose symmetry axis coincides with the bond 
vector joining monomers 13 and 12 with its tip located 
at monomer 12 having a tip angle ^max, (ii) the spherical 
shell centered on the same monomer 12 with inner and 



outer radii 



and 



This volume amounts to 



' '-'max 

)(r 



3 

max 



(5) 



The limiting values Tmin and 



are chosen such that 



U^'^irmin) = t/r'^lnnax) = SfceT and the extreme cone 
angle ^max is chosen such that ^'^''"'^(6'max) = SfceT. This 
implies that the volume Vc is the locus of points in space 
where the incremental increase in vibrational potential 
energy remains a few fceT. 

In the depolymerized state J, where the reacting 
monomer ii is shown in dark blue, it is still located 
close to the filament end monomer 12 but is now free 
(in the sense that it only interacts with the filament end 
through excluded volume interactions). More precisely, 
to be in the depolymerized reacting state J, the free re- 
acting monomer ii must lie in a spherical layer centered 
on monomer 12 with volume 



n3 ^ 



(6) 



where iimin is such that C/2(i?min) = SfceT and i?max = 
2 X 21/6 a (namely twice the range of the potential in 
Eq. ©). 

We now envisage the energy modifications involved by 
the transitions between / and J by focusing on the incre- 
mental energy of the reacting monomer ii with respect 
to the rest of the system. In state /, the total incre- 
mental energy, which contains both intramolecular and 
intermolecular contributions, can be written as 



tot,/ CV.7 
^.1 = % ' 



vib,/ 



£0, 



(7) 



where the contributions from the stretching energy 
Ui^^{r) and the bending energy ^t"^"^ involving the re- 
acting monomer have been regrouped into e^'^". The sum 
of all excluded volume pair interactions U2(r) between 
monomer ii and the rest of the system is written as e°^. 

The total incremental energy of monomer ii in state 
J is the sum of all its excluded volume pair interactions 
with the rest of the system (at a different location in 
comparison with state /) , giving 



tot, J CV,J 



(8) 



Our chemical step model is best expressed in terms 
of reactive monomers. During the course of the simula- 
tion, any free monomer in the vicinity of the active end 
of a filament, i.e. located within the spherical layer of 
size Vs centered on the end- monomer of the filament, will 
be coined as "free reactive monomer" meaning that it is 
susceptible to polymerize from a J-type depolymerized 
state. Note that the reactive end of a filament is active 
unless the filament has reached the maximum size of z 
monomers, in which case it can only depolymcrize. 



Similarly, the last monomer at the active end of a fil- 
ament containing a minimum of four monomers (trimers 
are not allowed to dissociate as filaments must have at 
least three monomers) will be considered as a reactive 
end-monomer, i.e. susceptible to depolymcrize from an 
/-type polymerized state, as long as it is located within 
the low intramolecular energy volume Vc- 

During the time window in which a particular 
monomer is reactive, it can be the object of a chemical 
step. The occurrence of such an event being sampled in 
a Poisson distribution of times with nominal frequencies 
fixed (for reasons which will be justified later) to 

^dcpoi ^ ^i^j ^ Gxp(-/3eo), (9a) 



^pol ^ ^J->7 ^ ^ 



Vs + Vc 
Vc 



(9b) 



The frequency prefactor v is a free parameter which can 
be used to tune the trial frequencies (to tune the effec- 
tive barrier height of the reaction) and hence the chemi- 
cal rates in the system without affecting the equilibrium 
state. 

When a chemical step involving a particular reactive 
monomer is selected by Poisson sampling statistics to oc- 
cur at a specific time t, depending on the nature of the 
reactive monomer, an MC attempted move (/ J) or 
( J ^ I) is made by sampling its new position uniformly 
from the volume Vs or Vc relative to the new J or / state. 
This new configuration is then accepted or rejected on the 
basis of an acceptance probability chosen to be 



cv, J 



P/-'^ = Min[l,exp(-/3(e^; 



P. 



ib.7 



n ))], (10b) 



^.1 

cv. J 

- e 



where (3 = (fceT)^^- In practice, an explicit sampling of 
this acceptance probability finally decides the success of 
the chemical step at the time t originally decided from 
the Poisson samplings of reaction times. When the chem- 
ical step is accepted, the reacting monomer is instanta- 
neously and definitively transferred to its new location 
and topological status while keeping its linear momen- 
tum unaltered. The dynamical trajectory is then con- 
tinued for times greater than t on the basis of the new 
microscopic state dynamical variables. On the contrary, 
if the attempted chemical step is not accepted, no reac- 
tion is recorded and the reactive monomer simply follows, 
with the rest of the system, the dynamical trajectory it 
would have followed in the absence of any chemical step 
at time t. 

Wc note that finally, whether a chemical step is finally 
recorded or not, the reactive monomer remains reactive 
as long as it has not diffused out of the relevant Vc or Vs 
reactive volume in which it lies at times just later than t. 
This point illustrates the coupling between monomer dif- 
fusion and reactive events. It indicates that the effective 
kinetic rates will not be strictly proportional to v at high 
frequencies as a too large value of the attempt frequencies 
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will necessary lead to ineffective sequences of polymeriza- 
tion/depolymerization steps involving the same filament 
end and reactive monomer. 

The justification of the above Poisson/MC algorithm 
to deal with (de)polymerization steps requires the proof 
that it satisfies detailed balance at equilibrium for any 
pair of specific microscopic states {I, J). This can be seen 
easily by considering the ratio of the number of transi- 
tions between the two states iV^~^'^ and N'^'^^ which can 
be expressed as 

NJ^^ exp(-/36*f"^):..^->^(l/K)P/er^ 

where both the numerator and the denominator contain 
four factors equal or proportional respectively to (i) the 
probability to be in the starting microscopic state / or 
J, (ii) the probability per unit time to make an attempt 
starting from the original state towards the alternative 
state of the reactive monomer, (iii) the probability to 
reach the specific final microscopic state J in 14 or / 
in Vc in the attempted move and finally (iv) the prob- 
ability to accept this final state, hence to record a suc- 
cessful chemical step. To verify that the above ratio is 
unity, one simply needs to envisage explicitly two cases, 
whether the combination e^^'^ + ej''^'^ — e^^'"^ is positive 
or negative. Making it explicit for the former case, use 
of Eqs. 0, dH]), dH) and Cni) leads to 

N'^' ^ oxp - e^f '"]) exp(-/3eo) ^ ^ 

N'^^' exp(-^[e-'Vef'^-e->'^]) 

where e*°*'^ (or similarly for J) denotes the total incre- 
mental energy of the reacting monomer in state /. We 
note that alternative schemes are possible. However we 
observed that the present one maximizes the acceptance 
probabilities and thus optimizes the efficiency of the al- 
gorithm. 

The MD algorithm is a step by step procedure where 
h is the MD time step and At (usually chosen such that 
At = Ush where Us is an integer) is the time interval be- 
tween two successive stochastic collisions in the MPCD 
procedure. The chemical reactions are integrated to the 
hybrid MD/MPCD scheme in the following way. At any 
discrete MD time t, we make an exhaustive list of reac- 
tive monomers. This consists of filament-end monomers 
as well as free monomers which may appear more than 
once (in the list) if they are active with respect to more 
than one filament end. Then the occurence of each re- 
active step (/ — > J) or {J — )■ /) during the next time 
interval between t and t + ft, is independently sampled by 
taking a random number ( from a uniform distribution 
between and 1 and compared to 1 — exp {—v^^^'h) or to 
1 — exp [—v'^^h). Any detected occurence is recorded. 
In most cases, after scanning the full list of potential re- 
action steps, we find zero or one occurence: in the latter 
case, the chemical step algorithm is applied at time t until 
final MC acceptance or rejection of the new location and 



new topological link of the reactive monomer. In prac- 
tice, we have never observed more than two occurrences 
during a MD step h and in the rare double occurence 
cases, we attempt the reactions in succession. 



D. Simulation ensemble including chemical reactions 

Theoretical developments and simulation experiments 
will be performed in the (TVt, TV/, V, T) ensemble in which, 
besides the volume {V) and the temperature (T), the to- 
tal number of monomers Nt and the total number of fila- 
ments Nf are both independently fixed. This ensemble is 
subjected to the additional constraint that the filament 
size, expressed in number of monomers, is restricted be- 
tween a minimum of three and a maximum of z. To 
any macroscopic state specified by a set of fixed ther- 
modynamic variables, there will correspond an equilib- 
rium free monomer number density pi and a filament 
length distribution {pi}i=3,z which depends on the set 
of equilibrium constants Ki associated to the relevant 
(de)polymerization steps connecting filaments of length 
i and z -I- 1 . Dynamic fluctuations at equilibrium involve 
various diffusive and reactive aspects which are in prin- 
ciple also species dependent. Considerable simplification 
occurs if the equilibrium constant K can be treated as 
filament size independent. The length distributions then 
become simple exponentials and the concept of (an over- 
all) average polymerization rate, U , per active filament- 
end and the corresponding average depolymerization rate 
W per active filament-end can be employed. In a sub- 
critical regime {U < W) or in the supercritical regime 
{U > W), the length distribution must be respectively 
a decreasing or an increasing function of the filament's 
size, as a result of detailed balance. The ratio U /W is 
equivalent to the ratio between the actual free monomer 
density pi and the state point dependent reference ("crit- 
ical") density, pc = 1/K. This ratio pi = pi/pc = U/W 
directly indicates whether the state point is subcritical 
(pi < 1) or supercritical (pi > 1). Within an ideal so- 
lution treatment, an i independent equilibrium constant 
turns out to be a function of only the temperature, say 
K^(T), and a concentration independent critical density 
p° = \/K^ more directly separates subcritical (pi < p°) 
and a supercritical (pi > pi) concentration regions. 

We close this section with some considerations about 
our choice to impose maximal and minimal bounds to 
the filaments in the context of living biofilaments, sub- 
ject to single monomer polymerization or depolymeriza- 
tion. The presence of a lower bound to the size of ex- 
isting filaments and the absence of consideration of any 
explicit nucleation mechanism imply that the number of 
filaments Nf is strictly constant and fixed by initial con- 
ditions. As biofilaments are often generated by initiator 
proteins (like profilin for F-actin^) in vivo and in vitro 
biomimetic experiments, the imposition of a lower bound 
size in our simulations can be interpreted as a way to link 
the number of filaments to a fixed number of permanently 
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active initiators. The imposition of an upper bound to 
the filament length follows primarily from our wish to 
treat living biofilaments in a simulation box for a broad 
range of equilibrium conditions, both in subcritical and 
supercritical conditions. In the former case, the upper- 
limit will have marginal effects if z is large with respect 
to the characteristic size associated to the decaying expo- 
nential distribution. An equilibrated situation in super- 
critical conditions cannot be maintained with free living 
filaments as they would grow indefinitely. A stationary 
distribution is only possible if some confinement effect 
stops the unlimited growth of filaments. The z upper 
limit should thus be seen as a simplified way to impose a 
confinement constraint allowing the establishment of the 
equilibrium state in supercritical conditions, a situation 
after all symmetric of the role of the lower bound (three 
monomers) in subcritical conditions. We finally note that 
in simulations, it will be pragmatic to adjust the size of 
the box L such that the maximum length of the filaments 
satisfies z/q < L. 



III. THEORETICAL FRAMEWORK 

We provide below some theoretical predictions on the 
filament length distribution and on the dynamics of fil- 
ament length fluctuations for the equilibrium system 
at fixed temperature described in the previous section. 
Our system consists of living filaments with size re- 
stricted between lower and higher bounds, subject to 
single monomer (de)polymerization steps at imposed to- 
tal monomer number density and fixed global filament 
number density. The theoretical framework provides the 
tools for a consistency check of our hybrid MD simulation 
technique by comparing in the next section the results of 
simulation experiments at various state points with the 
theoretical predictions. We start by approaching the fil- 
ament distribution at a macroscopic level but rely to sta- 
tistical mechanics to discuss the microscopic expressions 
of the equilibrium constants. The second part of this sec- 
tion deals with the time evolution of the filament length 
populations on the basis of kinetic equations adapted to 
our system and we provide the microscopic expressions 
of these chemical rates. 



A. The filament length distribution and equilibrium 
constants 



Restricting ourselves to equilibrium, we denote the free 
monomer density by p\ and the filament densities by 
ipi)i=3,z- The constraints on the system are 

Pi + 3/93 + 4/94 + ... + zp, - pf = 0, (13) 



and 



P3 + P4 + ••• + p. - p/ = 0, 



(14) 



where pt and pf are respectively the fixed monomer and 
filament global densities. 

Chemical equilibrium requires that for each 
(de)polymerization reaction Ai + Ai ^ ^i+i (im- 
plying filaments of length i and i + I and a monomer), 
the chemical potentials of the different species satisfy 
fii+i — fj-i — fj.1 =0 which can be transformed into an 
expression relating densities 



piPi-i 



(15) 



where the equilibrium constant Ki is a state point de- 
pendent quantity. Under the reasonable assumption that 
these equilibrium constants beyond the trimer (i = 3) can 
be treated as independent of the filament size, the indi- 
vidual equilibrium densities can be expressed explicitly 
in terms of a unique constant K which then turns out to 
carry the whole state point dependence of the distribu- 
tion. 

Applying pi = pipi^iK recursively and combining 
with Eqs. (|13|) and (jl4p . one gets an implicit expression 
linking the monomer density pi and K, 



Pt = Pi + Pf 



ELAKpi) 



i-3 



EUiKpi) 



i-3 



(16) 



and an expression for filament densities in terms of pi 
and 



Pi = Pf^ 



(KpiT 



(17) 



Analytical expressions corresponding to Eqs. (|16|) and 
(fT7)) can be obtained. Using reduced densities pt = PtK, 
Pf ~ Pf^ ^^'^ Pi ~ Pi^ to simplify expressions, we 
exploit the properties of finite geometrical series Sn = 
ELo(Pi)' = (1 - (/3i)"+')(l - Pi)~' (note that Sn is fi- 
nite and continuous around pi ~ I for n non negative and 
finite) and r„ = ^{PiY — PidSn/dpi. Equations 

p6p and p7p can then be rewritten as 

Pt = Pi + Pf- 



S 



z-3 



Pi +Pf 



3-2pi-(/9i)^-^[l+z(l-/9i) 

(l-pi)(l-(/5i)^-2) 



(18) 



and 



, (1 -/3i)(pi)» 3 

Pf— TTTI^^ = C'exp(y (i -3)), (19) 



l-(pi) 



where y = lnpi and C = p/(l — pi)[l — {pi)^~^]~^. 

The average size of the filaments in terms of monomers 
is then given by 



[3-2pi-/3r'(l + ^[l-pi])] 

(i-/9i)(i-pr') 



(20) 



Let us stress that Eqs. ([T51) - P^ . shown respectively in 
Figs, d] and O can alternatively be established via a 
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FIG. 4. pi is plotted vs pt with p; — 0.6699 (fixed) after 
numerical evaluation and inversion of Eq. (jlSp . The solid line 
is for z — 18, the dashed line for z = 40 and the dash-dotted 
line for z = 72. It is seen that when 2 — ^ oo, pi — >■ 1 as pt 
becomes large. The green (light) line at pi = 1 is a guide 
for the eye. Such curves have only a universal character for 
ideal solution conditions as the reduction factor is only 

a function of temperature (see text). 



statistical mechanics route (see Appendix . A semi- 
grand canonical ensemble with similar constraints takes 
into account all possible arrangements of Nt monomers 
distributed into Nf filaments and A'^i remaining free 
monomers, at fixed volume and temperature. Such an ap- 
proach is richer as it allows the calculation of additional 
properties like the composition fluctuations in finite sys- 
tems. Let us further point out that Eqs. (fT5|) -(fro |l are not 
closed equations, unless K can be established separately. 
The reduction factor (l/K) being generally state point 
dependent, the above relations are not universal curves. 
They merely represent consistency relationships between 
pi, Pi and K, for an arbitrary state point specified by (T, 
P/> Pt)- 

The free monomer density pc = l/K which corre- 
sponds to pi = 1 is known as the critical density. If 
an explored state point in our ensemble accidently yields 
pi = Pc, all filament length densities are equal, as seen 
in Eq. IH]), and iVav = (3 -I- z)/2. As discussed in 
the next subsection on the dynamics of the filament 
length, the polymerization and depolymerization steps 
are then equiprobable and each filament performs a one- 
dimensional (ID) discrete random walk bounded from 
below at i = 3 and from above at i = z. 

For an arbitrary state point, pi will differ from unity. 
If pi > 1 (supercritical conditions) , the filaments are sub- 
ject to a similar but biased bounded random walk where 
polymerization steps are more frequent than depolymer- 
ization steps (except for i ~ z) and the stationary fila- 
ment length distribution is an increasing function of the 




FIG. 5. Theoretical bounded filaments distributions for sub- 
critical (pi < 1) and supercritical (pi > 1) monomer densities 
shown as P'"'(i) = pi/pf vs i. The solid, dotted, dashed and 
dash-dotted lines are for pi — 0.6,0.9, 1.11 and 1.67 respec- 
tively in the conditions of the plot for 2 = 18 in Fig. 3] The 
solid horizontal line corresponds to the pi = 1 case. 



length in order to satisfy detailed balance. The oppo- 
site situation is observed for subcritical conditions, as 
pi < 1. Figures|3]and[5]illustrate respectively the equilib- 
rium relationship Eq. (fT5|) and the distribution Eq. (|19p . 
Equation psp for finite upper bound z cannot easily be 
inverted and must be solved numerically. 

Provided pi < 1, the z ^ 00 limit of the upper value 
of the filament range can be taken in Eqs. ([IS]), ([TO]) 
and (|2(jp which simplify respectively to 



Pt Pi + Pf 



(3 - 2pi) 



p, =p/(l-pi)(pi)*-3 



with 



Here Eq. (|21|) can now be inverted to yield 
1 



(21) 
(22) 
(23) 



Pi = 



{l + pt-2pf) 



(l + pt-2p/)2-4(pt-3p/) 



(24) 



All of the above expressions starting with Eq. pB]) . in 
particular the exponential distribution of filament sizes 
in Eq. ()19p . apply to equilibrium situations where the 
various chemical (de)polymerization steps are character- 
ized by a unique (size independent) equilibrium constant 
K. Such a simplified macroscopic approach will be found 
to be fully relevant to our model system. Therefore, we 
now exploit standard statistical mechanics expressions of 
the equilibrium constants, first within an ideal gas ap- 
proach and then within a dilute solution approach where 
interaction effects, between filaments and monomers, are 
included at the second virial coefficient level. 
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1. Ideal solution 

We first consider an ideal solution case, reminiscent 
of the reacting ideal gas treated in standard booksf^i^ 
where all interactions between filaments and free 
monomers are neglected, but where (de)polymerization 
steps nevertheless take place between the different (phan- 
tom) species, the filament + monomer solute system be- 
ing coupled to a heat bath "solvent" at fixed tempera- 
ture T. In that case, the equilibrium number densities 
(pi, {pi}i=3,z) must satisfy the chemical equilibrium re- 
lationship 



PiPi-i 



(25) 



where K^, function of the temperature only, is related 
to single filament or single monomer canonical partition 
functions by- - 



AT = V 



(26) 



The right hand side (RHS) of Eq. ^ can now be 
derived analytically for our filament model Hamilto- 
nian provided only vibrational intramolecular interac- 
tions Ui'^^{r) and (f>°™'^{9) arc considered. As discussed 
in Section III Al such a situation is valid as long as 
Lc — Ip > Iq, because the intramolecular "long-range" 
monomer-monomer interactions which in fact act only 
at short physical distances cannot contribute for slightly 
bending rods. Using Eqs. ([T]) and (jl]) for i > 2, the 
kinetic contributions cancel out and thanks to the inde- 
pendence of the individual vibrational contributions, one 
finds an i independent equilibrium constant for the RHS 
of Eq. (gg) to be 



A'"(r) = 27rexp(/3eo) / dO sin 6* cxp [-/30'^™'^(6')] 



X / dr r^exp[-(3U^'^{r) 
Jo 



(27) 



27rexp(/3eo) 



[1 - exp(-2/3fc')] / lo 



(1 + erf [xo]) (1 + xl) + xo exp {-xl/2) 

(28) 



Here xq = {Phf^lo, /3 = [kBT)-^ and 

erf[s] = y/2fn ( dt cxp(-t^). 

JO 



(29) 



The derivation of Eq. (^5)) from Eq. (P7)) is exact and 
relatively straightforward, requiring the use the integral 



/ duw^exp(-M^/2) = 

J a 



1 b 



— erf(u) — u exp (— u^/2) 



In our simulations, the values of the intramolecular pa- 
rameters are chosen such that individual bond lengths 
have negligible fiuctuations [xq 3> 1) and filaments semi- 
flexible [Lc < Zp). Therefore, Eq. of the equilibrium 
constant simplifies to 

A'O(T) = (27r)3/2/3-3/2A:-i/2(^')-i;2 exp(/3eo) (31) 

which more directly shows the impact of the model pa- 
rameters on the equilibrium constant. 



2. Effects of interactions 

At this point, we introduce the activity coefficients /i 
and {/i}i=3,2 which measure for monomers and filaments 
of different lengths, the deviation from ideal behavior 
in the composition term of the chemical potential. In 
terms of these, the equilibrium constant associated to 
the reversible depolymerization of a filament of length i 
can be expressed as^ 



Pi 



flfi- 



Pi Pi 



±K\T) = 



(32) 



In Eq. p2[) . Ki is a function of the thermodynamic state 
(and no longer a function of only the temperature). In 
presence of intermolecular forces, the i dependence of 
Ki for filaments is expected to be very weak while the 
magnitude of the intermolecular effects on AT^, for our 
model filament model, should primarily be related to the 
change in covolumc resulting from the addition or the 
removal of a monomer at the end of a filament. This point 
can be analyzed quantitatively within the framework of 
the low density expansion of the activity coefficients. In 
Appendix |B1 we obtain a first-order estimate of Ki / A'° 
as 



fifi- 



- = 1- [2feii + - hii] pi - 

Z 

[bik + (1 + 5k,i-i)bk,i-i - (1 + 5k,i)bki] Pk 

(33) 



fc=3 



(30) 



Oip') 



where the b factors are the expansion coefficients of the 
pressure in powers of the activities of the various species 
in the mixture. These factors can be expressed as follows 



(34a) 
(34b) 
(34c) 
(34d) 



where the pairs of indices imply corresponding two body 
interactions and the pair configurational integrals Z^v 



bii = 


— {Zii 




bli = 






bu = 


-{Zu - 


v\ 
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between species u and v arc defined as 



Zuu — 



Quv ^ 

{€)' 



when u ^ 



(35) 
(36) 



In these expressions, 5° and are one body canonical 
partition functions and Quv and Quu are two body canon- 
ical partition functions for different and similar bodies 
respectively. 

To be specific in this calculation of the second virial 
coefficients, we treat our filaments as rigid rods and we 
concentrate on the expressions relevant to the first order 
term in pi in the Ki/K^ expression given by Eq. (|33p . 
We have^ 



Z„ = V^y drexp(-/3un(r)), (37) 
Zu = V J drexp(-/3ui,;(r,0,a;)), (38) 

where un is the monomer-monomer pair potential and 
uii is the intcrmolccular potential between a filament of 
length i with center of mass at the origin and with (arbi- 
trary) fixed orientation a; and a free monomer located at 
r. The explicit expressions of the corresponding b factors 
are 

foil = ^ y" dr [cxp i-l3un{r)) - 1] , (39a) 
bu = [ dr [cxp {~/3uu{r, 0, w)) - 1] . (39b) 



These integrals can be estimated by approximating the 
free monomers as hard spheres of diameter a and fil- 
aments as hard sphcrocylindcrs built as a cylinder of 
height H = {i — 1)Iq and radius (t/2, terminated by two 
hemispheres of radius a/2. Performing the integrals gives 



foil = -vrv2 
foi. = -vr 



(40) 



= -[- + (z-l)]^a2/o, (41) 



where V^"'^ and V^^°'^ are the excluded volume of respec- 
tively a sphere or a spherocylinder, when approached by 
another sphere. Supposing that in Eq. (|33p the monomer 
density is dominant with respect to the other terms in- 
volving filament densities, we get the net covolume effect 
on the equilibrium constant as. 



K/Ko - 1 



" 2, 



(42) 



This result gives at least the order of magnitude of inter- 
molecular effects on the equilibrium constant and more- 
over suggests, as expected, that the i dependent effects 
should be small. 



B. Mean field rate equations and rate constants 

1. Effective (de) polymerization rate constant 

The (de)polymerization processes which take 
within the system imply the following reaction 



fcoff 



place 



(43) 



where Ai represents a free monomer, Ai a filament con- 
sisting of i monomers and Ai-^-i one with i + 1 monomers. 
In the mean field kinetics model, the i independent fcon 
and fcoff are respectively the polymerization and the de- 
polymerization rate constants for filaments undergoing 
reactions at their active ends. Chemical equilibrium im- 
plies that the mean field equilibrium constant K , a func- 
tion of the thermodynamic state, satisfies K = fcon/^off- 
Let us first consider polymerization events in an equi- 
librium situation with densities (pi, {pi}i=3.z). The com- 
bination konPiPi expresses the number of polymerization 
steps per unit of time and per unit volume involving 
specifically filaments of length i. In our general statis- 
tical mechanics model the filaments can only fluctuate 
between i — 3 and i = z. Accordingly, the mean of the 
total number of reactive pairs for polymerization per unit 
volume (at any time) is given by V* pi{p3 + p4 + ... + pz-i) 
where the volume factor equals to 



V* 



An 



dr r^gomir) 



(44) 



in terms of the equilibrium radial pair distribution gom{r) 
between any filament active end and surrounding free 
monomers. In the ideal solution case, it simplifies to 
V* = Vs- According to our microscopic kinetic model, 
the total number of successful polymerization steps per 
unit volume and per unit time can be written as a prod- 
uct of three contributions, namely, the total number of 
reactive free monomers per unit volume, the polymeriza- 
tion attempt frequency Eq. (j9bp . and finally, the poly- 
merization acceptance probability Eq. (jlObp . We thus 
have 



pol 



--konPi{pz + Pi + ■■ + Pz-l) = 
V*pi[p3 + Pi + ■■ + Pz-l) X I' 



Vc. 



V. + V. 



X (Min[l,exp(-/3((e: 



cv,/ 



vib.7 



^cv. J\\l\co 



)W 



(45) 



Here the average (...)'^° runs over all polymerizing events 
J I sampled uniformly in the reacting volume V^- 
The incremental excluded volume energy of the reacting 
monomer n in the different states / or J are represented 



Also ( 



■ib.7 



is the non-negative incremen- 



tal vibrational energy of the reacting monomer n in state 
/ obtained by grouping the stretching potential Ul^^ir) 
and the bending potential 0'^°"'^(6') provided by Eqs.([T|) 
and (g]). 
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Now considering the depolymerization events at equi- 
librium, the average number of reactive (breaking) pairs 
per unit volume is given by f*{p4 + P5+P6 + --- + Pz) 
where the factor /* is the fraction of end monomers at 
the active end of the filaments which lie in the reactive 
volume Vc- 

The total number of successful depolarization per unit 
of time and per unit of volume is, using depolymeriza- 
tion frequency Eq. (|9a[) . and depolymerization accep- 
tance probability Eq. (jlOap . given by 



of equations governing population dynamics reads 



^dopol ^ + P5 + Pe + -Pz) 

{P4 + P5 + P6 + --Pz)!* X V- 



■ exp(-/?eo) 



Vs + Vc 

X (Min[l, cxp (-/3(e-'^ - {^'^ + e^^'^'^Mr. 

(46) 

In this case, the average {...y^ runs over depolymerizing 
events I ^ J sampled uniformly in the spherical layer of 
volume Vs. 

To summarize, we finally get the rate expressions 

fcon = 



X (Min[l,exp {-^7'' + ^T^') " C'''))])^°, 

(47) 



■cxp(-/3eo) 



X (Min[l,cxp (-/3(6-.-' - (e-'^ + e^'^-'MY^- 

(48) 

The order of magnitude of these rates for reasonably di- 
lute solutions can be estimated by neglecting excluded 
volume interactions, using /* ~ 1, ^ Vs and the ideal 
solution equilibrium constant A'" given in Eq. (|3ip . as 



h 



off 
I cstim 



z^exp(-/3eo), 
A'°i^cxp(-/3eo)- 



(49a) 
(49b) 



dP(3,t) 

dt 
<lP{k,t) 



AP{z,t) 
dt 



~?7P(3,i) + VKP(4,t) (50a) 

-{U + W)P{k, t) + UP{k - 1, t) 
\-WP{k + l,t) (50b) 

'WP{z,t) + UP{z-l,t) (50c) 



where P{i,t) is the probability to observe filaments with 
length i at time t and in Eq. (j50bp the index k runs from 
4 up to z — 1 . This set of equations is solved for constant 
rates U and W for given initial conditions P(i,0). This 
approach which treats independent filaments is physically 
justified if the monomer density (apparently uncoupled 
to the dynamical scheme) is assumed to be homogeneous 
in space and in time. Well stirred systems with com- 
pensating external addition or removal of free monomers 
provide plausible non-equilibrium stationary conditions 
for which the set of Eqs. (|50p is pertinent. 

In our context we explore different specific values of 
U ~ fconPi and W — kos. The Eqs. (|5D|) can be tested by 
our simulations if we probe the dynamical filament length 
fluctuations at equilibrium. For example, by picking up 
a subset of filament of length j at some time t = and 
following the time evolution of their size distribution, one 
should get an evolution from the original delta function 
back to the equilibrium distribution. This distribution 
is simply the conditional probability P{j, t; i) that a fil- 
ament having a length i at an initial time t = has a 
length j at time t. 

This conditional probability time evolution will reflect 
the statistical properties of the bounded biased random 
walk induced by the chemical steps. However, if the 
starting size i is not immediately close to a boundary 
at i = 3 or at i = z, the behavior of P{j,t;i) will be 
representative of an unbounded biased random walk as 
long as the measuring time t remains much smaller than 
the (shortest) mean first passage time towards one of the 
boundaries. At long times, boundaries have a major in- 
fluence (in any subcritical or supercritical regime) on the 
re-establishment of the equilibrium distribution. 



2. Mean-field filament length dynamics 

The filament length dynamical evolution is usually 
treated theoretically in terms of a chemical master equa- 
tion with the reaction rates treated as input constants. 
For biological filaments like actin, filaments change by 
one unit and, if we disregard for the present time 
the existence of variants in the complex formation of 
the monomers and the modification of ATP into ADP 
through hydrolysis, only the polymerization rate U and 
the depolymerization rate W are relevant for the single 
monomeric species considered. In our case where the fil- 
ament lengths fiuctuate between i = 3 and i = z, the set 



IV. ADOPTED MICROSCOPIC PARAMETERS AND 
LIST OF SIMULATION EXPERIMENTS. 

All the physical quantities will be expressed in a system 
of units based on: the length ag of the side of the cell 
in MPCD, the energy kBT (where fee is the Boltzmann 
constant) and the mass rus of a solvent particle. The 
resulting time unit is expressed as tq = aQy/mTITk^T)- 

For the monomer/filament model, we adopt a unique 
set of parameters inspired by previous works on semi- 
flexible chains in MPCD solvent^^i^^ or in BD solvent^!, 
namely beads of mass m = 5.0 with size a = 0.44545 
and energy parameter e = 3. Within the fllament. 
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TABLE I. Table with simulation results for different values of A''t. The particle packing fraction is represented by rj. We fix 
Nf — 80 in all the experiments with their size limited between 3 and z — 18. Here pf™ is the reduced monomer density obtained 
from the slope of the logarithm of the equilibrium filament length distribution, and /of™ is the measured single monomer density. 
For comparison, we also tabulate pi numerically calculated from Eq. (|18|) . The polymerization (U) and depolymerization (W) 
rates are calculated directly by counting. 



Exp. 


Nt 


V 


Pi 


„sim 

Pi 


Pi 


U 


= feoff 


fcoii ^ U/pl 


Gl 


1076 


0.0121 


0.60658 ± 0.004 


0.12203 ± 0.0002 


0.12227 


0.00244 


0.00399 


0.02000 


G2 


1744 


0.0196 


0.91490 ± 0.003 


0.17994 ± 0.0006 


0.18009 


0.00356 


0.00389 


0.01978 


G3 


2008 


0.0225 


1.00756 ± 0.001 


0.19805 ± 0.0003 


0.19808 


0.00388 


0.00385 


0.01959 


G4 


2338 


0.0262 


1.13775 ± 0.001 


0.22141 ± 0.0003 


0.22168 


0.00432 


0.00379 


0.01951 


G5 


3311 


0.0372 


1.78321 ± 0.008 


0.33824 ± 0.0002 


0.33831 


0.00637 


0.00360 


0.01883 


G6 


4000 


0.0449 


2.39206 ± 0.020 


0.44891 ± 0.0001 


0.47836 


0.00821 


0.00345 


0.01829 



monomers are linked by hard springs with equilibrium 
length Iq = 0.5 and force constant k = 1600. The bend- 
ing parameter is set to k' = 20, which implies a persis- 
tence length of k' /{k-sT) ~ Ip/lo = 20. These values lead 
to a reactive volume Vg = 3.81855 (see Eq. ^) defined 
by i?min = 0.44545 (for which the potential is SksT), 
and Rmax = 1, while those defining the reactive volume 
Vc = 2.9 X 10-2 (see Eq. ©) are given by r^in = 0.43876, 
r^^^ = 0.56124 and = 0.55481 rad. 

The MPCD parameters are chosen to match a well 
studied state point ^2 given by a solvent density of psoiv = 
5, box collision angle of 130° and a collision frequency 
1/At = 10. At this solvent state point, it is known that 
a particle of mass m = 5ms interacting at infinite dilution 
via MPCD collisions with solvent, has a diffusion coeffi- 
cient^^ of Dq = 0.043 so that an estimate of the diffusion 
time is tq = (0.463428)Vi:»o = 5.0 (where r = 0.463428 
is the distance where the potential U2{r) = k^T). 

With the adopted kinetic parameters v = 5.0 and eo = 
6.9 one gets at ksT = 1, the ideal solvent equilibrium 
constant K'^ = 4.88373. This allows the estimates of 
rates, from Eqs. (glD, as fcj;|'™ 5. x 10"^ and /c^^*™ = 
2.5 X 10"^ respectively. 

All the simulation experiments, conducted at a com- 
mon temperature /cbT = 1, follow Nf ~ 80 filaments 
with size fluctuating between three and 2 = 18 which 
are enclosed in a cubic volume V = 18^ containing 
29,160 MPCD solvent particles. Note that the box 
length (L = 18) is twice the maximum allowed filament 
length. Periodic boundary conditions are applied in all 
three dimensions. Table [T] reports the total number Nt of 
monomers which is specific to each experiment, together 
with experiment codes Gl - G6 increasing with Nt and 
the associated volume packing fraction 77. In the sim- 
ulations, the MD integration time step h, associated to 
the Verlet algorithm, has been set to ft- = 0.002, which 
corresponds exactly to At/50. 

A typical simulation setup involves an equilibration 
run of 5 X 10^ corresponding to « 2.5 x 10'^ depolymer- 
ization steps per active filament end or to « IO^tu. Four 
runs of time duration of 10^ are then performed in suc- 
cession to get four independent estimates of the averages 
of interest which are then used to determine the global 




FIG. 6. Plot of P'"^{i) = pi/pf vs i for two different pi 
values, namely subcritical experiment G2 (black circles) and 
supercritical experiment G4 (black squares). The solid black 
lines are the corresponding fits. The two dashed straight lines 
indicate the expected distributions in ideal solution conditions 
for the same state point (Nt, Nf, V, T). 



average (over a total time of 



4. X 10^) and an 



estimate of the associated statistical error. 



V. SIMULATION RESULTS 

A. Equilibrium distributions 

The experiments differing in the total number of 
monomers Nt in the box lead to a resulting free monomer 
density (see Table H]) and a distribution of filament 
densities for all accessible lengths (monomer number). In 
all cases, the distribution of filament lengths turns out to 
be exponential, confirming the hypothesis of a unique 
(state point dependent) equilibrium constant K for all 
filament lengths. In Fig. [SI P°'^{i) = pi/pf (in log scale) 
versus i (for two cases) is plotted. It shows a linear be- 
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FIG. 7. The equilibrium constant calculated from the sim- 
ulations. The open black circles are calculated from the 
free monomer density and the exponential filaments densi- 
ties while the red stars are estimated from the ratio of the 
rates {1/ p{)U /W . Within the error bars, the agreement be- 
tween the two estimates is excellent. The solid red line is the 
ratio (l/pi)/(pi)/<7(pi) estimated from the fits. The dotted 
line represents the equilibrium constant for the ideal case. 

havior with a slope a providing directly the reduced free 
monomer density pi = piK ~ exp(a). All the pi data 
obtained from the filament distributions of the various 
experiments are gathered in Table [H It shows that as Nt 
increases, the system goes from subcritical to supercriti- 
cal conditions, with the experiment G3 being very close 
to the critical situation. In Fig. [HI we also show the dis- 
tribution expected (dashed lines) for the same state point 
{Nt, Nf, V, T) if ideal solution conditions were applicable 
with equilibrium constant computed from Eq. (PT|) . 

By dividing pi by pi for each line of the Table HI one 
gets a first estimate of the equilibrium constant K which 
is shown in Fig. [3 as a function of pi itself. We postpone 
the discussion of the increase of K versus pi until we 
confirm these data with a second estimate of K extracted 
from the number of chemical (de)polymerization steps in 
the next subsection. 



B. State point dependence of rates and related 
equilibrium constant 

The effective rates W = kos and U ~ konPi can be esti- 
mated, on the basis of the assumption that they are inde- 
pendent of filament length, using expressions in Eqs. (|45|) 
and (|46p . The quantities 7iP°' and n'^°P°^ in these expres- 
sions are obtained by the counting of the total number 
of successful polymerization and depolymerization steps 
observed during the simulation, divided by the volume 
V and the total simulation time. The rates for the six 
experiments are gathered in the Table |T] and are also re- 



q|^_, 1 , 1 , 1 , 1 , 1 

0.1 0.2 0.3 0.4 0.5 
Pi 

FIG. 8. Polymerization rate U (circles) and depolymerization 
rate W (squares) vs pi for all experiments. The solid lines 
are linear (W) and quadratic (U) fitting curves discussed in 
the text. In the inset, the global growth rate J — U — W is 
indicated (diamonds). The solid line showing J = is only 
for reference. 

ported in the Fig. [51 

We observe that W slightly decreases with increasing 
pi, a property which must be related to a decrease in 
the rate of acceptance of an attempted depolymerization 
step due to higher packing density. Similarly, the poly- 
merization rate is not strictly linear in pi and a negative 
quadratic term can be extracted from a fit. We perform 
the three parameters combined fit W{pi) = f{x) = ax+b 
and U{pi) = koapi = g{x) = cx^ + K'^bx as we im- 
pose that for pi — >■ 0, the ratio of fcon and k^ff repro- 
duce the thermodynamic value of if ° = 4.88373. We find 
a = -0.00165663, b = 0.00417753 and c = -0.00461497. 
To first order, it gives the relationship K/K^ = 1 + Ax 
with A = c{bK"y^ - ab'^ = 0.169. This value can 
be compared with the theoretical prediction given by 
Eq. (glD, which gives A' = 0.112 if the size of the 
monomer is taken to be equal to the distance at which 
the energy is of the order of unity {a w 0.463). The es- 
timate has the correct order of magnitude but is about 
35% too low perhaps because the covolume effects of the 
filaments are included in the fitted quantities {U and W) 
but are absent in the theory. We believe however that 
the main source of increase of the equilibrium constant 
with pi finds its origin in the covolume effect described 
in Eq. dUl). 



C. Dynamic fluctuations of filament length 

At equilibrium, we consider the (static) population 
probability P'^'^{i = 3, z) and the conditional probabil- 
ity P{j, t; i) for a filament to have a length j at time t if 
its length was i at time i = 0. In terms of these, one can 
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FIG. 9. C{t) vs t for experiments Gl (squares), G3 (circles) 
and G5 (diamonds). The solid curves represent the theoretical 
expressions of C{t) based on Eqs. (|51|l and (|52p where the 
conditional probabilities are computed by solving according to 
standard methods, the mean field population dynamics given 
by Eqs. (fSU)) . using the U and W actual tabulated values of 
the same state point. 




200 



FIG. 10. (AA'^(f))jv(o)=i vs t at short times for different 
values of pi. The lines refer (from bottom) to an average over 
filaments i = 8, 9 and 10 in experiment Gl, an average over 
filaments i = 9, 10 and 11 in experiment G3 and an average 
over filaments i = 11, 12 and 13 in experiment G5. 



quantity measured and shown in Fig. [TO] is 



express the time correlation function of filament length 
Nhy 



(iv(o)iv(i)) = EE*j'^''W^(j'^;*)- 



(51) 



=3 3=3 



The normalized correlation function of the length fluctu- 
ations around their mean is formally 



Cit) 



{N{Q)N{t)) - {Ny 



(52) 



where {N) = iVav and {N^) = ELs Figure [9] 

shows the variation in C(t) with t for the experiments 
Gl, G3 and G5. These curves are in excellent agree- 
ment with the theoretical expressions of C{t) based on 
Eqs. ([5T|) and ([5^ where the conditional probabilities 
are computed by solving the mean field population dy- 
namics Eqs. (|50p . using the tabulated values of U and W 
at the same state point. We observe that the relaxation 
time of these fluctuations is maximum at critical density 
(experiment G3), as a result of the uniform distribution 
of filament lengths which requires a filament to explore a 
wide range of lengths before loosing memory of its initial 
value. 

In Fig.[TOl we observe the growth or decay of the length 
of filaments in various super or subcritical regimes, focus- 
ing on the short time behavior of the conditional proba- 
bilities for particular filament lengths i which are present 
in reasonable number and located relatively far from the 
boundaries at three or z at the same time zero. The 



(AiV(t))^(0)^, = ^ jP(j,t;z) 



(53) 



i=3 



and is also averaged over a few appropriate filament 
lengths. Since the largest (de)polymerization rates are 
of the order of 0.01, it can be expected that for times 
^obs = 100 ~ 200, the filament populations of P{j, t; i) 
for J = 3 or j = z are still negligible and the behavior of 
{AN{t)) is representative of a situation at the same den- 
sity pi of unbounded filaments. We have verified that the 
slopes of the straight lines observed in Fig. [TO] are in full 
agreement with the global growth or decay rate J shown 
in the inset of Fig. [5] 



VI. DISCUSSION 

The methodology presented in this work applies to the 
sampling of a reactive (semi-grand) canonical ensemble 
consisting of a mixture of free monomers and a fixed num- 
ber of self-assembled linear polymeric filaments. Specif- 
ically, the filaments undergo reversible single monomer 
end-filament polymerization/dcpolymerizations at one or 
both active ends of any filament. Regarding the filament 
model, a discretized wormlike chain is used as the basis 
of an extension towards a "living wormlikc chain" ver- 
sion where the contour length is variable. In this section, 
we mention some specific details of our method, some of 
them with respect to alternative approaches proposed in 
the literature. We briefly indicate how our methodology 
can be generalized to include additional or alternative 
kinetic processes found in biofilamcnts such as filament 
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rupture, distinction between ADP and ATP complexes 
and irreversible hydrolysis processes. We organize this 
discussion by regrouping some purely methodological as- 
pects first, and then consider the explicit case of actin to 
illustrate the versatility of our approach. 



A. Methodological aspects 

The dynamic trajectory of our solute (reactive mixture 
of filaments and free monomers), coupled to a bath of sol- 
vent particles, is obtained via a numerical scheme which 
incorporates MC reactive steps in the hybrid MD-MPCD 
method. For non-reacting fluids, the hybrid algorithm 
conserves the total energy but in our case, the MC chem- 
ical steps give an overall canonical ensemble character to 
the sampling. The imposed canonical temperature y^^" 
is set in the acceptance probabilities of the chemical step 
Eqs. (jlUp and in the attempt depolymerization frequency 
Eq. ([5a|) ). Between two successive chemical events, the 
system's total energy is conserved but it is discontinu- 
ous when a chemical step is successful. Pragmatically, 
the initial equilibration run to probe a new state point 
is started with initial velocities for solute monomers and 
solvent particles such that the average kinetic energy per 
degree of freedom is set equal to As the sys- 

tem equilibrates while undergoing chemical steps at the 
same prescribed temperature, the instantaneous average 
kinetic energy per degree of freedom fluctuates in time 
around the expected stable value. The choice of cou- 
pling of the new chemical steps proposed in this paper to 
a MPCD solvent may be questioned. This (continuous 
space) particulate description of the solvent locally con- 
serves momentum and hence preserves long wavelength 
fluid hydrodynamics as well as allows the boundary con- 
ditions to be easily adjusted. This may play an important 
role in highly conflned geometries such as in biofilaments 
in a microchannel or actin bundle developments in fllopo- 
dia. 

The stochastic scheme to model the 
(de)polymerization steps is a purely instantaneous 
local perturbation. This is meant to represent at a 
mesoscopic level (coarse graining in space and time), 
the direct coupling effect of many internal and external 
degrees of freedom of the reacting pair "fllament end"- 
"protein complex" which are eliminated. The modeling 
of the chemical step introduces a local discontinuity 
in the position of the reacting monomer but not in 
the monomer velocity. This is the price to be paid in 
order to succeed in devising a scheme in which the free 
reacting particle can be effectively moved to/from a 
position in space where it participates with the other 
degrees of the filament to an equilibrated intramolecular 
configuration. In real systems, as one protein binds 
to an existing protein filament in solution, it requires 
some time to modify its internal structure, to adjust 
its orientation and its translational location so that 
intermolecular interactions between units stabilize the 



newcomer. This waiting time corresponds in our model 
to the sampled reaction time in the Poisson process 
governed by the attempt frequency in Eqs. ((HJ, which 
can be tuned via the free parameter v. The reactive steps 
are modeled as local events obeying micro-reversibility. 
They are effectively coupled to the solute degrees of 
freedom through structural properties like the filament 
end- monomer/free- monomer pair correlation function, 
and through the various filament and free monomer 
diffusive processes. Our simulation results indicate a 
weak variation of the rates with state point at moderate 
solute packing fractions when the effective kinetic rate 
constants are computed for a fixed set of parameters 
regarding the stochastic model (reacting volumes Vc and 
Vs, and attempt frequency v). We have pointed out that 
the order of magnitude of the rates can be estimated 
from the ideal solution behavior (Eqs. P^)) ). 

There is a large amount of literature on simulation 
methods/studies of flexible or semiflexible equilibrium 
polymers but most of them resort to various purely MC 
schemes which sample a grand-canonical ensemble (with 
addition/removal of individual molecules homogeneously 
in the simulation box) subject to additional constraints 
on chemical potentials in order to impose chemical equi- 
librium^!^ or use sometimes artificial MD schemes to 
compute only the equilibrium static quantities 4i Our ex- 
ploration of the reactive semi-grand canonical ensemble 
by MD gives similar information on thermodynamic and 
structural quantities (such as distributions of filament 
length populations, various monomer/filament pair cor- 
relation functions etc), but in addition yields dynamic 
information on diffusive and reactive processes both at 
equilibrium, as illustrated in the present work, and po- 
tentially outside equilibrium. 

Reaction-diffusion aspects of biofilaments (actin) in 
presence of explicit reactive free monomers have been 
studied recently in inspiring BD studiesi^"— The 
depolymerization steps are modeled as explicit end- 
monomer detachments directly triggered by an indepen- 
dent stochastic process with fixed depolymerizing rate. 
The polymerization step relies on the detection of free 
monomers entering, by diffusion, into a capture zone 
volume defined at the reactive end. Subsequent to its 
entry, the monomer is automatically treated as reactive 
and subsequently attached to the pre-existing flexible fil- 
ament. In all these studies, the opposite polymeriza- 
tion and depolymerization steps do not satisfy micro- 
reversibility. While Lee and Li u^^'^^ studies are very 
specific and concerned with far from equilibrium net- 
work developments close to a moving disk (including 
branching and capping phenomena) to investigate the 
origin of disk motility, Guo et al. are concerned with 
more basic phenomena on single F-actin filaments re- 
garding size and composition fluctuations. They treat 
a single unbounded fllament surrounded by free Brown- 
ian monomers at a strictly constant concentration, mim- 
icking an homogeneous stationary non-equilibrium sit- 
uation. Purely ADP-actin complexes have been con- 
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sidered in the first study concerned with a stationary 
growth or shrinkage of the filament. In the next work, 
treadmilling has been produced by a similar BD study 
in presence of three types of actin complexes namely, 
ADP-, ADP+Pi- and ATP-actin. (Here Pi represents 
inorganic phosphate group). In these BD studies, for 
pragmatic reasons, chemical reactions are artificially ac- 
celerated with respect to a free monomer diffusive time 
scale of reference, by increasing both kinetic rates and 
free monomer concentration in order to follow diffusion 
and chemical events within a unique simulation time win- 
dow allowing significant statistics on the slowest dynam- 
ical events. Our methodology is different but works at 
the same mesoscopic level and therefore, if applied to 
the actin case, it would face the same wide spectrum of 
time scales problem discussed above. Progress on fur- 
ther time-space rescaling needs to be envisaged such as 
considering the elementary monomeric unit treated in 
the model to represent a collection of a few ATP-protein 
and/or ADP-protein complexes, with an effective slower 
diffusion constant and adapted chemical steps. 

Coming back to the large classes of equilibrium and 
nonequilibrium self-assembled polymers or filaments con- 
sidered in the introduction, the kinetic processes which 
are involved in their self-assembly and disassembly are 
known to be of very different nature. Cylindrical micelles 
can break anywhere and reform by end-recombinations 
and this mechanism is potentially plausible to some ex- 
tent in biofilaments also4^ Biofilaments are usually sub- 
ject to ATP-hydrolysis while the building blocks of the 
filaments are ATP-protein complexes which can thus be 
transformed in other complexes when integrated into fila- 
ments. In our method, we only restrict the chemical steps 
to be single unit end polymerization or depolymerization 
and that too to single monomer species. Our approach 
can however be easily extended to cover all these cases 
as illustrated in the next section where we discuss some 
aspects of actin modeling. 



B. Parametrization of our generic biofilament model to 
F-actin 

For illustrative purposes, we discuss below the appli- 
cation of our method to actin at the level of one free 
monomer particle being one single protein complex with 
ADP or ATP. Referring to the study of actin by Guo et 
q1^^7j28 discuss in detail the parametrization needed 
in the equilibrium polymer case where any monomer, 
either free or integrated in a F-actin filament, remains 
permanently in the actin-ADP complex state. We then 
briefly discuss extension towards non-equilibrium actin 
filaments subject to irreversible ATP hydrolysis. 

To begin with, we remark that the diffusion time of 
free monomers in our study, depending on the mass of 
the free monomer and various parameters of the MPCD 
solvent, turns out to be 5 tq. This can be used to set the 
basic time unit as tq = 0.5 x 10~^ s when the experimen- 



tal Globular-actin diffusion coefficient of 10~^^ m^/s is 
taken into account4^ The size of the monomer is 0.5 qq, 
fixing our length unit to oq = 10 nm. As actin has a 
double-stranded helix structure;^ it is usually in simpli- 
fied models as a single strand such that each monomer 
addition increases the effective contour length of the fila- 
ment by half the size of a uniti^ This implies we should 
take Iq = 0.25 ag. The force constant k is not a sensible 
physical parameter but it should be fixed to a sufficiently 
large value to avoid unphysical fluctuations of the contour 
length. The bond fluctuations are small with respect to 
In if the quantity a;o defined in Eq. ((^ satisfies ^ 1. 
In the present discussion, we can keep the value k = 1600 
in simulation units, as this value leads to xq = 10. Next, 
we must fix k' such that we reproduce the experimental 
value of the persistence length of the pure ADP-F-actin 
at T = 293 K which is = 9 /im.— This gives a value of 
k' = 3.6 X 10^ kBT. The value of the equilibr ium constant 
will next be imposed in the model by using Eq. (pij) 
valid as long as ^ 1, and by exploiting the freedom 
in choosing the bonding energy parameter eg. We thus 
find Co = 18.7 kBT in order to get the experimental crit- 
ical concentration of pc = = 1.8 /iM4^ These set 
of parameters still require the rate to be adjusted to 
match the experimental data of the off-rates (the on-rates 
are automatically satisfied if K'^ is correctly fixed). As 
one has different values at both ends, namely k^g = 5.4 
s^^ at the barbed end and k^^ = 0.25 s^^ at the pointed 
end)^ one gets by applying Eq. (|49a|) . i^*^^ w 7 x 10^ s~^ 
and i^P° « 3 X 10^ s'^. 

Moving now to F-actin filaments consisting of three dif- 
ferent forms of complexes, namely ATP-actin, ADP+Pj- 
actin and ADP-actin, requires associating a flag to each 
monomer in the model. Any active end of the filament 
("barbed end" or "pointed end") may differ in monomer 
composition (hydrolysis status) and hence specific rate 
constants for each case must be adjusted in the model 
parameters. The kinetic model for the associated ir- 
reversible hydrolysis of the ATP form needs also to be 
specified, as well as the associated kinetic rates. To dis- 
cuss the parametrization in this context, let us adopt 
a simplified treatment^^ where the intermediate form 
of ADP-|-Pi-actin is merged with the ADP-actin form 
so that only two distinct forms are retained and where 
the hydrolysis step, which irreversibly changes the ATP- 
complex in the filament into an ADP-complex, takes 
place only at the boundary of an ATP cap and pure ADP 
section (vectorial process). In such a block-copolymer, 
the different parameters fc'^^^,e^°^ and i^^^^ for the 
ADP block at the pointed end (and at the barbed end 
in case the ATP cap has vanished) can be adjusted on 
those of the pure ADP case discussed previously. One 
can proceed in the same way for the ATP block adjust- 
ing the parameter k'^"^^ to match its stiffer character 
(persistence length of 15-17 ^rr4^i^) and further impose 
values of e^"^^ , v^"^^ for the barbed end (and the pointed 
end in case the ADP block would have vanished) in or- 
der to match the corresponding experimental values of 
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kinetic rates and hence of the equilibrium constants. Fi- 
nally, an independent stochastic process with rate fixed 
to its experimental value would convert the first ATP- 
complex of the cap into the ADP form. Actin simula- 
tions would be run as explained in Section |TT] but with a 
much large variety of reactive monomers including ATP 
(or ADP) barbed end (or pointed end) monomers, ATP 
monomers at the diblock boundary and ATP (or ADP) 
free monomers in the reactive volume centered on an ac- 
tive end (barbed or pointed). For the hydrolysis reaction, 
if sampled during a MD step, the reaction is automati- 
cally performed as an irreversible step. In all other cases 
where a monomer is selected to perform a specific reac- 
tion, an attempted MC move is performed which is then 
subject to the appropriate acceptance criterion. 



such that the total number Nt of monomers is fixed. It 
reads 



Q{Nt,Nf,V,T) 



^ iVi! 

A'i,(Afi)i=3,.;{C} 



(Al) 

where qi is the canonical partition function of a single 
filament of size i or of a free monomer for i = 1 ,^1^ Ow- 
ing to the fluctuating size of the filaments between their 
lower and upper limits, the sum in the partition function 
involves all possible arrangements of Nt monomers which 
satisfy the explicit constraints, jointly denoted formally 
as {£} in Eq. ([JT]) . 



iVi + 3A^3 + 4Ar4 + ... + zN^ - iVt = 0, 



(A2) 



VII. CONCLUSION 

To conclude, we have presented a hybrid MC-MD 
model to simulate living biofilaments in the framework 
of the widely used wormlike chain model, with reactive 
dynamics at the ends giving rise to contour length vari- 
ations. We believe that our mesoscopic approach, be- 
ing formulated within a statistical mechanics framework, 
could be the starting point for numerous applications and 
extensions. The study of the growth of filament bundles 
in free or in confined space are in progress. More gener- 
ally, in the multi-scale problem of living biofilaments dy- 
namics, our methodology offers a well defined model at 
an intermediate level which should open natural bridges 
to more atomistic models on one side and on the other 
side, to more coarse-grained models to be based on fur- 
ther systematic and controlled rescaling procedures. On 
a side note, our simulation approach could also be easily 
adapted to study polymer synthesis in good solvent by 
living polymerization, in particular, to follow the time 
evolution of the size distribution of the polymers 1^ 
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Appendix A: Statistical mechanics treatment 

The partition function of our semi-grand canonical re- 
active ensemble treats independent solute entities divided 
into a set of N f living filaments with polymer size fluc- 
tuating through single monomer (de)polymerization be- 
tween ? = 3 and i = z and a set of A^i free monomers 



and 



(A3) 



In addition, we suppose that the state point equilib- 
rium constant K associated to monomer attachment 
to/detachment from a filament is independent of the fil- 
ament size. 

This theoretical treatment can be justified within two 
approaches which are both relevant to the interpretation 
of our simulation experiments involving interacting and 
reactive filaments: 

• It strictly applies to the ideal solution case (sub- 
jected to constraints given by Eqs. (|A2[) and (jA3|) ) 
as the equilibrium constant K'^{T), independent of 
the filament size for our model, depends only on 
the temperature (see Eq. ([5T|) )). 

• It provides an approximate treatment based on a 
mean-field approach whereby the partition function 
is written in terms of independent effective single 
filament partition functions qi. The qi and hence 
the related equilibrium constant K supposed to be i 
independent, are state point dependent. The state 
point dependence of K needs to be approximated 
to close the set of density equations (see Eq. (01])). 

The recursive application of Eq. with Ki — K gives 
be rewritten as 



so that the partition function can 



Q(iV„iV„F,T) = ,(--^-V^(| 



E 



{Nt-3Nj) 



1 



Afi,(JVi).=3, 



Ni\N3\...NJ \ V 



K 



-Ni 



(A4) 



The average number densities of monomers and filaments 
in the thermodynamic limit can be estimated^E by search- 
ing for the largest term of the partition function Q. This 
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requires solving the constrained global minimum 



In the above equations, the series are limited to quadratic 
terms in Zi. These relations can be inverted by writing 



dNi 



In 



d 



1 (K 

1 



dN, 



■In 



- A = 0, (A5a) 
(i = 3,z), (A5b) 



zi 



Pi + aiiPi + ^ aiipipi 



(B5) 



i=3 



where A and p are Lagrange multipliers related respec- 
tively to the constraints Eqs. (|A2[) and (|A3[) . Use of the 
Stirling's approximation leads to the set of equations 



ln7Vi+ln(^) + A = 0, 


(A6) 


ing 




In Ni +iX- fi = [i 


= 3,z). (A7) 


an 


= -26ii, 






an 




In terms of reduced number densities 


p — pK, one gets 


a-ii 




pi = cxp (-A) 


(A8) 




— —bki- 




Ni = exp {-[iX + p\) = p\ exp (-/i) 


(i = 3,z). (A9) 


Using Eq. (|B2b|, one g 


ets 



Zi = pi + cii^pf + aiipipi + ^ aktPiPk (B6) 

where the coefficients a^- can be obtained by substitu- 
tion of the latter expansions in the density expansions 
(Eq. (jB4[) ). followed by term by term identification, giv- 



(B7) 
(B8) 
(B9) 
(BIO) 



The combination of Eqs. (jA3p and (|A9p gives exp (— /i) = 

Nf (X]i^=3 Pi) I hence leading to the filament densities 
Eq. (jl7p . while the implicit equation in free monomer 
density Eq. (fTBl) follows by direct substitution of filament 
densities in the constraint Eq. (|A2[) . 



/i = 1 - 2biipi - biiPi, 

i=3 



/j = 1 - "^hiPi - bupi - bkipk, 

k—3^k^i 



(Bll) 

(B12) 



Appendix B: First-order estimate of equilibrium constant in 
the presence of excluded volume interactions 

The pressure expansion of a mixture of monomers and 
filaments in terms of activities reads^^ 



zi+^Zi + biizl + ^ bi,z1 + ^ buZiZi 

i—3 i— 3 i— 3 

z — 1 z 

+ E E b,,^^^l+0{z^), (Bl) 



(B13) 

The substitution of these density expansions of the ac- 
tivity coefficients in the equilibrium constant expression 
(Eq. ([211)) leads to the requested expression (Eq. ([33])). 



where the activities are defined by 

zi = f,p,=q",cxpiPp,)/V, (B2a) 
z,; = /,p, =q°exp(/3^0/^^- (B2b) 

Here the coefficients bij are related to two-body integrals 
by Eqs. and Eqs. The number densities can 

then be expressed as 



Pi = Zl 



Pk = Zk 



/dpp\ ^ ^^(^ ^ 26nzi + ^H^O, (B3) 



1=3 



V ^-^1 / 

i-Q^j^Zk{l + 2bkkZk + bikZk+ E ^''i^j)- 

(B4) 
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